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Motivated by a recent application of the Kolmogorov-Johnson-Mehl-Avrami (KJMA) model to 
the study of DNA replication, we consider the one-dimensional version of this model. We generalize 
previous work to the case where the nucleation rate is an arbitrary function I(t) and obtain analytical 
results for the time-dependent distributions of various quantities (such as the island distribution). 
We also present improved computer simulation algorithms to study the ID KJMA model. The 
analytical results and simulations are in excellent agreement. 
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I. INTRODUCTION 

Consider a tray of water that at time t = is put 
into a freezer. A short while later, the water is all 
frozen. One may thus ask, "What fraction f(t) of wa- 
ter is frozen at time t > 0?" In the 1930s, several scien- 
tists independently derived a stochastic model that could 
predict the form of f(t), which experimentally is a sig- 
moidal curve. The "Kolmogorov-Johnson-Mehl-Avrami" 
(KJMA) model [1 has since been widely used by 
metallurgists and other materials scientists to analyze 
phase transition kinetics 0. In addition, the model has 
been applied to a wide range of other problems, from 
crystallization kinetics of lipids 5], polymers 0, the 
analysis of depositions in surface science 0, to ecolog- 
ical systems and even to cosmology 0. For further 
examples, applications, and the history of the theory, see 
the reviews by Evans 01 > Fanfoni and Tomellini , and 
Ramos et al. [IT|. 

In the KJMA model, freezing kinetics result from three 
simultaneous processes: 1) nucleation of solid domains 
("islands"); 2) growth of existing islands; and 3) coales- 
cence, which occurs when two expanding islands merge. 
In the simplest form of KJMA, islands nucleate anywhere 
in the liquid areas ("holes"), with equal probability for 
all spatial locations ("homogeneous nucleation"). Once 
an island has been nucleated, it grows out as a sphere 
at constant velocity v. (The assumption of constant v is 
usually a good one as long as temperature is held con- 
stant, but real shapes are far from spherical. In water, 
for example, the islands are snowflakes; in general, the 
shape is a mixture of dendritic and faceted forms. The ef- 
fect of island shape - not relevant to the one-dimensional 
version of KJMA studied here - is discussed extensively 
m 0.) When two islands impinge, growth ceases at the 
point of contact, while continuing elsewhere. KJMA used 
elementary methods, reviewed below, to calculate quan- 
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FIG. 1: Definitions. In the KJMA model, a hole is the liq- 
uid domain between the growing solid domains (island). The 
island-to- island is defined as the distance between the centers 
of two adjacent islands. 



tities such as f(t). Later researchers have revisited and 
refined KJMA's methods to take into account various ef- 
fects, such as finite system size and inhomogeneities in 
growth and nucleation rates 0, . 

Although most of the applications of the KJMA model 
have been to the study of phase transformations in three- 
dimensional systems, similar ideas have been applied to a 
wide range of one-dimensional problems, such as Renyi's 
car-parking problem [T^ and the coarsening of long par- 
allel droplets [13j . Recently, we have shown that the one- 
dimensional KJMA model can also be used to describe 
DNA replication in higher organisms Briefly, in 

higher organisms (eukaryotes), DNA replication is initi- 
ated at multiple origins throughout the genome. A repli- 
cated domain then grows symmetrically with velocity v 
away from the replication origin. Domains that impinge 
coalesce. And finally, each base in the genome is repli- 
cated only once per cell cycle. Thus, if one views repli- 
cated regions as "solid," unrcplicated ones as "liquid," 
and the initiation of replication origins as "nucleation," 
all of the essential ingredients of the KJMA model are 
present. The purpose of the present two papers, then, is 
as follows: Here, in paper I, we discuss how to generalize 
the KJMA model for biological application. In particu- 
lar, we consider the problem of arbitrarily varying origin 
initiation rate (equivalent to arbitrarily varying nucle- 
ation rate in freezing processes). Then, in paper II, we 
discuss a number of subtle but generic issues that arise in 
the application of the KJMA model to DNA replication. 
The most important of these is that the method of anal- 
ysis runs backward from the usual one. Normally, one 
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FIG. 2: Kolmogorov's method, (a) Spacetime diagram. In 
the small square box, the probability of nucleation is IoAxAt, 
where 7o is the nucleation rate. In order for the point X to 
remain uncovered by islands, there should be no nucleation 
in the shaded triangle in spacetime. (b) Kinetic curve for 
constant nucleation rate Io'- f(t) = 1 — exp(—Iovt 2 ). 



starts from a known nucleation rate (determined by tem- 
perature, mostly) and tries to deduce properties of the 
crystallization kinetics. In the biological experiments, 
the reverse is required: from measurements of statistics 
associated with replication, one wants to deduce the initi- 
ation rate I(t). This problem, along with others relating 
to inevitable experimental limitations, merits separate 
consideration. 

In the mid-1980s, Sekimoto showed that the analysis of 
the KJMA model could be pushed much further if growth 
occurs in only one spatial dimension 16]. Sekimoto used 
methods from non-equilibrium statistical physics to de- 
scribe the detailed statistics of domain sizes and spacings, 
as defined in Fig. ^ I n particular, he studied the time 
evolution of domain statistics by solving Fokker- Planck- 
type equations for island and hole distributions, for con- 
stant nucleation rate I(i)=const. His approach has since 
been revisited by others (e.g. 0). 

Below, we extend Sekimoto's approach to the case of 
an arbitrary nucleation rate I(t). As mentioned above, 
this case is relevant to the kinetics of DNA replication in 
eukaryotes. We also present two algorithms to simulate 
ID nucleation and growth processes that are both much 
faster than more standard lattice methods 1181. 



II. THEORY 

A. Island fraction f(t) 

We begin with the calculation of f(t), the fraction of 
islands at time t in a one-dimensional system. We write 
as f(t) = 1—S(t), where S(t) is the fraction of the system 
uncovered by islands (i.e., the hole fraction). In other 
words, S(t) is the probability for an arbitrary point X at 
time t to remain uncovered. If we view the evolution via 
a two-dimensional spacetime diagram [Fig. [2(a)], we can 
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which has a sigmoidal shape, as mentioned above [see 

Fig. EDO] ■ 

We note that Kolmogorov's method can be straight- 
forwardly applied to any spatial dimension D for arbi- 
trary time- and space-dependent nucleation rates I(x, t). 
Similar "time-cone" methods can yield f(t) in t he p res- 
ence of complications such as finite system sizes [lj, [l^ . 
Unfortunately, this simple method cannot be used to cal- 
culate the distributions defined in Fig. ^ except that it 
can partly help solve the time-evolution equation for the 
hole-size distribution (see below). 



B. Hole-size distribution ph(x,t) 

We define ph(x,t) as the density of holes of size x at 
time t. For a spatially homogeneous nucleation function 
I(t), the density ph will also be spatially homogeneous 
(The hole size x should not be confused with the genome 
spatial coordinate X). The time evolution ph{x,t) then 
obeys 

dp h (x,t) dp h (x,t) 

= 2v - I(t)xph{x,t) 



dt 



dx 

/>oo 

-2/(t) / p h (y,t)dy, 



(3) 



where v is the growth velocity of islands and lit) is the 
spatially homogeneous nucleation rate at time 1 116| . The 
first term on the right-hand side describes the effects on 
Phix,t) of domain growth in the absence of coalescence 
and nucleation. The second term accounts for the anni- 
hilation of a hole of size x by nucleation, while the last 
term represents the splitting of a hole larger than x by 
nucleation. Eq.OJwas solved by Sekimoto for J(i)=const., 
while Ben-Nairn et al. derived a formal solution for ar- 
bitrary I(t) [2l]]. Below, we show that the solution of 
Bcn-Naim et al. can also be obtained directly by apply- 
ing Kolmogorov's argument. In Fig. [3| we see a hole of 
size x flanked by two islands. In order for such holes to 
exist at time t, there should be no nucleation within the 
parallelogram ABCD in the spacetime diagram. Similar 
to the calculation of the hole fraction S(t), we obtain the 
"no nucleation" probability in the parallelogram as 
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FIG. 3: Spacetime diagram. The hole-size distribution 
ph(x,t) is proportional to the probability po(x,t) for no nu- 
cleation event occurs in the shaded parallelogram ABCD (see 
text). 



where g(t) = J Q I(t')dt'. The domain density n(t) and 
the hole fraction S(t) are related by definition as follows: 

(5) 
(6) 



n W = / Ph(x,t)dx 
Jo 



S(t) — / xph{x 1 t)dx. 
Jo 



Since the hole-size distribution ph(x,t) is proportional 



to po(x,t), we can write ph(x,t) = c(t) ■ po(x,t). By 
integrating this equation and using Eq. [5] we obtain 
c(t) = n(t) ■ g{t)/S{t). Putting this back into Eq. El 
we obtain an equation for n(t): 



1 dn(t) 
n(i) dt 



-2v.g(t) + *-§-. 



(7) 



This is a first-order linear equation and can be solved 
exactly. Using the boundary condition n(0) — 1, we solve 
Eqs.|3and|3to find 



n(t) 
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(t) ■ e ~ 2v H 9(t')dt' . 



Ph ( X ,t) = g(t) 2 -e- 9 ^ X - 2v ^9(t')dt\ 



(8) 
(9) 



These are just exponential functions of x, with decay 
constants that monotonically decrease as a function of 
time. 



C. Island distribution pi(x,t) 

In analogy to Eq. [3] and following [lfj , the time evo- 
lution of the island distribution pi{x,t) is governed by 
drift, creation, and annihilation terms, as follows: 



dpi(x,t) 
dt 



-2v 



dpi(x,t) 
dx 



I(t)S(t)6{x) + 2v 



Ph(0,t) 
n(t) 2 



pi(x - y, t)pi(y, t)dy - 2n{t)p i {x 1 1) 



(10) 



Again, the first term on the right-hand side represents 
the effects of domain growth. The second term accounts 
for the creation of islands of zero size, initially. [S(x) is 
the Dirac delta function.] The last two terms represent 
the creation and annihilation of islands by coalescence, 
respectively. We note that the prefactor 2vph(0,t)n(t)~ 2 
can be obtained by writing it as a(t), applying dx to 
Eas. 151 and 1101 and then comparing the two. 

Unfortunately, we cannot solve Eq. 1 101 using the simple 
arguments that worked for ph{x,t). The main difference 
is that a hole is created by nucleation only, while an island 
of nonzero size is created by growth and/or the coales- 
cence of two or more islands. Thus, Pi{x, t) is given by an 
infinite series of probabilities for an island to contain one 
seed, two seeds, three seeds, and so on. Nevertheless, we 
can still obtain the asymptotic behavior of pi(x, t) for ar- 
bitrary I(t) by Laplace transforming the above evolution 
equation, as in [lfj. 

Applying J Q dxe~ sx to Eq. ^| we find 



dpi(s,t) 
dt 



-2v[s + 2g(t)]pi{s,t) (11) 



where f>i(s,t) = J Q °° e sx pi(x,t)dx, with initiation con- 
ditions pi(s, 0) = 0. We can further simplify Ea. ITTlbv 
defining Gi(s, t) = exp \2v f g(t')dt'] -pi(s, t), which then 
obeys 



dGj(s,t) 
dt 



-2v[s + g(t)]Gi(s, t) + 2vGi(s, tf + I(t). 



If we write Gi(s,t) as 



Gi(s,t) = s + g(t)+X(s, t), 



(12) 



(13) 



we find that X(s, t) obeys the (nonlinear) Bernoulli equa- 
tion El: 



dX(s,t) 

at 



= [ S + g(t)}X(s,t)+X( S ,t) 2 . (14) 



Solving Eq. 1141 and substituting back into Eq.^J we find 
the Laplace transform p~i(s, t): 
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(15) 



We cannot perform the inverse Laplace transform 
of the above equation, even for the simple case of 
J(t)=const. [i.e., g(t) ~ t] [Tol Il7j. However, from the 
form of denominator in Eg. 1151 we observe that pi(s,t) 
has a single simple pole along the negative real-axis at 
\s = s*(t)\ -C 1 for t » 1, regardless of the form 
that g(t) may have. Since the inverse Laplace transform 
can be written formally as the Bromwich integral in the 
complexjslane (i.e., as the sum of residues of the inte- 
grand 23]), a standard strategy for obtaining the asymp- 
totic expression of pi(x,t) for x 3> 1 is to expand p~i(s,t) 
around s*(t) (|s*(t)| <C 1) to lowest order. Following 
Sekimoto's approach, we define K(s, t) to be the denom- 
inator in Eq. 1151 which becomes 



pi(s,t) 



s + g(t) 



1 dK(s,t) 1 



2v dt K(s,t) 



Around s = s*(t), Eg. 1151 can be approximated as 
e -Io9(t')dt' dK(s*(t),t) 1 



Pi(s,t) 



-2v dt 

e -J*g(t')dt' ds *^ 



dK(s*(t),t) 
ds 
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dt s-s*(t) 
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From Eq. 1161 we arrive at the following asymptotic ex- 
pression for pi(x,t): 
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for x, t 3> 1. Now, both the prefactor and the exponent 
[the pole s*(t)] can be obtained very easily by simple 
numerical methods. On the other hand, an approximate 
expression for s*(t) itself can be found by first expanding 
K(s, t) in powers of st and then solving iteratively using 
Newton's method |2J|. The result is 



1 / J x 4Jf - Jo J 2 



2J 4 



(18) 



where 



J n = / e-loa(t')dt' T n dT _ 
Jo 

As we shall show below, Eq. El describes the behavior of 
Pi(x,t) accurately for x > 2vt. 




FIG. 4: Constraint plane S : (ii + «2)/2 + h = x. 



D. Island-to-island distribution pi2i{x,t) 

While most studies of ID nucleation-growth have fo- 
cused on ph(x, t) and Pi{x, t) exclusively, the distribution 
of the distances between two centers of adjacent islands 
[the island-to- island distribution Pi%i(x, <)] has important 
applications. For instance, whether homogeneous nucle- 
ation is a valid assumption cannot be known a priori. In- 
deed, in the recent DNA replication experiment that mo- 
tivated this work, the "nucleation" sites for DNA replica- 
tion along the genome were found to be not distributed 
randomly, a result that has important biological implica- 
tions for cell-cycle regulation |25| . 

In the ID KJMA model, Sekimoto has shown that a 
constant nucleation function Iq cannot produce correla- 
tions between domain sizes [16|. We speculate that the 
same holds true for any local nucleation function I(x, t), 
a conclusion that is also supported by computer simula- 
tion l2Sj . Assuming a local nucleation function, wc 
can write the formal expression for Pi2i(x,t) directly in 
terms of pi(x,t) and ph(x,t): 



Pi2i(x,t) = C 



{ii,h,i 2 }(iS 



Pi(ii,t)ph(h,t)pi(i 2 ,t)dS, (19) 



where S designates the constraint plane shown in Fig. 0] 
[S : (i± + i 2 ) /2-\-h=x\. The normalization coefficient c 
can be obtained easily from the relation J Q pi2i{x, t)dx = 
J °° pi(x,t)dx — J °° ph{x,t)dx — n(t). From Eo. 1191 and 
Fig.QJ it is easy to see that J °° pi2i(x, t)dx = c[n(t)] 3 , and 
therefore c = [n(i)]~ 2 . Since the full solution for pi(x,t) 
is not known, we cannot integrate Eq. 1191 However, we 
can still obtain an asymptotic expression for Pi2i(x,t) 
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using Eas. 151 and 1171 For x 3> 1, taking into account the 
constraint S, we find 

p i2i (x,t) ~ f e -\s'(t)\-ii-9(t)-h-\a-(t)\-i 2dS 
J {ii ,h,i 2 }&S 

~ e- g W x + e - 2 l s *«l*[- 1 + g (t)x - 2\s*(t)\x] . 

(20) 

As we shall show later, the bottom Eq. [201 is an excel- 
lent approximation for all range of x and time t. Note 
that the first term on the right-hand side has the same 
asymptotic behavior as the hole-size distribution ph{x, t), 
while the exponential factor in the second term comes 
from the product of island-size distributions ~ e~' s wr* 1 
and ~ e~' s (*)l' 12 . The asymptotic behavior of Pi2i{x,t) 
is dominated by ph(x,t) for / < 0.5 and by pi(x,t) for 
/ > 0.5 (see below). But at all times, we emphasize that 
Pi2i{x,t) is asymptotically exponential for large x. From 
the mathematical point of view, both pi(x, t) and ph(x, t) 
have exponential tails at large x, and the integral of the 
product of exponential functions again produces an ex- 
ponential. 

III. NUMERICAL SIMULATION 

Often, one has to deal with systems for which ana- 
lytical results are difficult, if not impossible, to obtain. 
For example, the finite size of the system may affect its 
kinetics significantly, or the variation of growth velocity 
at different regions and/or different times could be im- 
portant. In such cases, computer simulation is the most 
direct and practical approach. 

For one-dimensional KJMA processes, the most 
straightforward simulation method is to use an Ising- 
model-like lattice, where each lattice site is assigned ei- 
ther 1 or (or — 1, for the Ising model) representing is- 
land and hole, respectively. The natural lattice size is 
Air — vAt, with v the growth velocity. At each timestep 
At of the simulation, every lattice site is examined. If 0, 
the site can be nucleated by the standard Monte Carlo 
procedure, i.e., a random number is generated and com- 
pared with the nucleation probability I(t)AxAt. If the 
random number is larger than the nucleation probability, 
the lattice site switches from to 1. Once nucleation is 
done, the islands grow by Ax, namely, by one lattice size 
at each end. 

Although straightforward to implement, the lattice 
model is slow and uses more memory than necessary, 
as one stores information not only for the moving do- 
main boundaries but also for the bulk. Recently, Herrick 
et al. used a more efficient algorithm [13 . Specifically, 
they recorded the positions of moving island edges only. 
Naturally, the nucleation of an island creates two new, 
oppositely moving boundaries, while the coalescence of 
an island removes the colliding boundaries. 

For the present study, we have developed two other 
algorithms, which have improved both simulation and 
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FIG. 5: Comparison of simulation times for the three algo- 
rithms discussed in the text. Circles are used for the lattice- 
model algorithm, squares for the double-list algorithm, and 
triangles for the phantom-nuclei algorithm. For each system 
size, the number of Monte Carlo realizations ranges from 5- 
20, and the lines connect the average simulation times. The 
double-list algorithm is two to three orders of magnitude 
faster than the lattice algorithm, while the phantom-nuclei al- 
gorithm ranges from three to five orders of magnitude faster, 
depending on the number of time points at which one records 
data. The filled triangles show the fastest case, with only one 
time point requested, while the open triangles show the slow- 
est case, where data are recorded at each intermediate time 
step. 

analysis speeds by factors of up to 10 5 (Fig. |5J. The 
first of these, the "double-list" algorithm, extends the 
method of Herrick et al. 0|- The second of these, the 
"phantom-nuclei" algorithm, is inspired by the original 
work of Avrami Q. 



A. The Double-List Algorithm 

Fig. describes schematically the double-list algo- 
rithm. The basic idea is to maintain two separate lists of 
lengths: {i} for islands, {h} for holes |2£j. The second list 
{h} records the cumulative lengths of holes, while {i} lists 
the individual island sizes. Using cumulative hole lengths 
simplifies the nucleation routine dramatically. For in- 
stance, for times t ranging between r and r + At, the 
average number of new nucleations is N = I(r)AxAt. 
Since the nucleation process is Poissonian, we obtain the 
actual number of new nucleations N = p(N) from the 
Poisson distribution p. We then generate N random 
numbers between and the total hole size, namely, the 
largest cumulative length of holes h max (the last element 
of {h}). The list {h} is then updated by inserting the 
N generated numbers in their rank order. Accordingly, 
{i} is automatically updated by inserting zeros at the 



6 



(a) 



{i} 





{i} 


{h'} 





i 


h ' 


= h 


1 


ii 


hi' 


= h +h, 


2 


h 


h 2 ' 


= h +h[ +h 2 



{h'> -^V 




Before nucleation 





{i} 


{h'} 







h ' 


1 




h,' 


2 


J 2 


ha' 



After nucleation 



_ T 

[ !l Jj 

ho' h,' h 2 ' 




« 


{h'} 





io 


ho' 


1 


ii 


h," = x 


1 1 A 1 1 
Random number between and h 2 ' 


2 


i 2 =0 


h 2 ' = old h/ 


3 


i, = old i 2 


h 3 ' = old h,' 



become much more complicated because the individual 
hole sizes would have to be taken into account as weight- 
ing factors in distributing the nucleation positions along 
the template. 

Fig. [S] compares the running times for two different 
algorithms: the standard lattice model vs. the contin- 
uous double-list algorithm described above. We wrote 
and optimized both programs using the Igor Pro pro- 
gramming language [23, and they were run on a typ- 
ical desktop computer (Pentium P3, 900 MHz). For 
both, we used the same simulation conditions: timestep 
At = 0.1, nucleation rate J(i) = 10~ 5 t, and growth ve- 
locity v = 0.5. Note that the performance of the lattice 
algorithm is O(N), whereas the double-list algorithm is 
roughly N 1 - 5 ^ 2 for 10 5 < N < 10 7 . The main reason is 
that the double-list algorithm has to maintain dynamic 
lists {i} and {h}. This requires searching and remov- 
ing/inserting elements (as well as minor sorting), where 
each algorithm is linear, or roughly 0(N 2 ) in overall. 
However, the double-list algorithm performed almost 3 
orders of magnitudes faster than the lattice algorithm 
even at a system size of 10 7 , and we did not attempt to 
improve the efficiency further, for example, by using a 
binary search. 

Finally, the relative storage requirements for the lat- 
tice algorithm compared to the double-list algorithm can 
be estimated by the ratio Ni at /n max , where Ni at is the 
total number lattice sites per unit length and n max is 
the domain density. Equivalently, one may use l m i n / Ax, 
where l m i n is the minimum island-to-island distance and 
Ax the lattice size. Since one usually sets up the sim- 
ulation conditions such that l m i n 3> Ax, the double-list 
algorithm requires much less memory. 
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FIG. 6: Schematic description of the double-list algorithm, 
(a) Basic set-up for lists {i} and {h'}. Note that {h'} records 
cumulative lengths, (b) Nucleation. (c) Coalescence due to 
growth. 



corresponding places. If {h} were to record the actual 
domain sizes as {i} does, the nucleation routine would 



B. The Phantom-Nuclei Algorithm 

Fig. describes schematically the phantom-nuclei al- 
gorithm. The basic idea is to capitalize on the ability 
to specify when and where in the two-dimensional space- 
time plane lie potential initiation sites, in advance of the 
actual simulation. Thus, in Fig. Q the circles, which 
represent potential initiation sites, are laid down in a 
first part of the simulation. One then uses an algorithm, 
described below, to determine which of these potential 
sites actually initiates (these are denoted by open circles) 
and which cannot fire because the system has already 
been transformed (these "phantom nuclei" are denoted 
by closed circles). 

The principal advantage of the phantom-nuclei algo- 
rithm is that one can find the state of the system at a 
particular time t without having to calculate the system's 
state at intermediate time steps. If one is interested in 
only a small number of system states, then the method 
can be significantly faster than the double-list algorithm. 
The filled triangles in Fig.[S]illustrate a hundred-fold im- 
provement compared to the double-list algorithm (and a 
10 5 -fold improvement relative to the lattice algorithm). 
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FIG. 7: Schematic description of the phantom-nuclei algo- 
rithm. The figure shows the distribution of potential initi- 
ation sites in the spacetime plane. The open circles denote 
sites that do initiate, while the "phantom" filled circles, lying 
in the "shadow" of the open circles, do not initiate. 

On the other hand, if information is needed at every time 
step (or if the number of phantom nuclei is very large) , 
the algorithm slows. The open triangles in Fig. show 
a simulation where information is collected at each time 
step. The run time is comparable to the double-list al- 
gorithm. Because there is no sorting operation, a lin- 
ear time scaling is maintained. The phantom-nuclei and 
double-list algorithms thus cross over at about two mil- 
lion sites. 

The phantom- nuclei algorithm is performed as follows: 

1. One generates the potential initiation sites in the 
two-dimensional spacetime plane. At each time, 
this is done in the same way as in the double-list 
algorithm. The only difference is that here, the 
number of sites at any time is calculated over the 
whole length of the system (regardless of its state 
of transformation). One uses two vectors to store 
the position and initiation time of every potential 
site. 

2. One removes all initiation sites that are in positions 
that have already transformed before their initia- 
tion time. (They lie in the "shadows" in Fig. 0) 
Because the growth velocity is known at each time, 
it is straightforward to implement this. Briefly, one 
first sorts the potential initiation sites by space. 
Then for each potential site (indexed by i, with 
position Xi and nominal initiation time ti), one cal- 
culates the position of the right-hand boundary n 
at the reference time point t. This is given by 
Ti = Xi + v(t — ti) for each i. One then proceeds 
through the list r^. If < rj for any j < i, then 
site i is discarded. Finally, one repeats the analo- 
gous process moving to the left, using the left-hand 
boundaries £j = Xi — v(t — ti). 

3. One calculates the desired statistics at the reference 
time point. This time point is arbitrary. For the 
filled triangles of Fig. [S] it is the last time step of the 
transformation process, while in the open triangles, 



it occurs at the next time step of the double-list 
simulation. (For the latter case, the statistics were 
then repeatedly calculated at each time interval.) 

In conclusion, we note that both the double- list and the 
phantom-nuclei algorithms are significant improvements 
on the more straightforward lattice algorithm. For simple 
initiation schemes, where it is possible to give a function 
I(t) for the inflation sites, the phantom-nuclei algorithm 
will generally be preferable. For more complicated initi- 
ation schemes, where the initiation of sites is correlated 
with the activation of earlier sites, the double-list algo- 
rithm may be preferable. In the next section, we present 
simulation results based on the double-list algorithm. 



IV. COMPARISON BETWEEN THEORY AND 
SIMULATION 

In Fig. [HI we compare the various analytical results 
obtained in the previous sections with a Monte Carlo 
simulation. Shown are ph(x,t), pi(x,t), and Pi2i{x,i) for 
I(t) = 1CP 5 £ at three different time points: t—50, 75, and 
100. The system size is 10 7 and the growth rate is v = 
0.5. The chosen form of accelerating I(t), linear in time, 
is the simplest nontrivial nucleation scenario. It is also 
relevant to the description of DNA replication kinetics 
in Xenopus early embryos, where the I(t) extracted from 
experimental data has a bilinear form [19j. 

The agreement between simulation and analytical re- 
sults is excellent. In particular, we emphasize that the 
analytic curves in Fig. [8] are not a fit. Note that, for 
x 3> 1, all three distributions decay exponentially as pre- 
dicted by Ea.|g|f!7l and ED] [The p h (x,i) distributions 
are simple exponentials over the entire range of x.] 

One interesting feature of Pi(x, t) is the inflection point 
in the interval < x < 2vt, where pi(x,t) is slightly 
convex. Such behavior is even more dramatic when 
/(i)=const. |16(, and pi(x,t) is strongly convex. In other 
words, pi{x,t) increases as x approaches 2vt~ , but sud- 
denly drops discontinuously at x = 2vt, decaying expo- 
nentially. This peculiar behavior of pi(x, i) originates 
from the fact that any island larger than 2vt must have 
resulted from the merger of smaller islands. Therefore, 
for x < 2vt, pi(x,t) has an extra contribution from is- 
lands that contain only a single seed in them, which 
makes pi(x,t) deviate from a simple exponential. Al- 
though such discontinuities are expected at every x — 
n ■ 2vt (n=l, 2, 3, . . .), higher-order deviations decrease 
geometrically and thus are almost invisible. 

Finally, the island-to-island distribution Pi2i(x,t) pro- 
vides important insight about the "seed distribution" 
and about the spatial homogeneity of the nucleation. 
Note that Pi2i{x,t) is not monotonic and has a peak 
at x > [see Fig. |S{c)]. This is not surprising be- 
cause Pi2i(x,t) — ► as x — > from Eq. ED On the 
other hand, we see that Pi2i(x, t) decays exponentially at 
large predicted in the previous section (Eq. I2U|) . In 
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(a) (b) (c) 

FIG. 8: (Color online). Theory and simulation results for I(t) ~ t. Size distributions are calculated at these timepoints: t = 
50, 75, and 100. (a) Hole-size distribution ph(x,t). (b) Island distribution pi(x,t). The inset plots f(t) vs. t, with the dot at 
t = 75 (/ = 0.5). (c) Island-to- island distribution pi2i(x,t). The analytical curves have been obtained by Eg. 1201 There is a 
crossover of the decay constant slightly after t=75 (/ = 0.5) (see text). The inset shows ph(x,t) and pai{x,t) for t = 50. All 
figures have the same vertical ranee of 10" 7 - 10" 3 (log-scale). 
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FIG. 9: (Color online). Decays constants of ph{x,t), pi(x,i), 
and Pi2i{x,t). The symbols are simulations, and the solid 
lines are theory. 



the corrections to this relationship in Eq. 20 imply that 
this holds true only for large x and t. Note that the actual 
minimum of Tin is at / > 0.5 because pi2i depends on pj 
and not pi alone. 

One final note about the island-to- island distribution is 
that, unlike Pi{x, t), it is a continuous function of x. The 
reason for this is that for any island-to-island distance x, 
the discontinuous pi(y < x,t) contributes to Pi2i(x,t) in 
a cumulative way, as can be seen in Eq. I19I This implies 
that there is no specific length scale where discontinu- 
ity can come in. From a mathematical point of view, 
this is equivalent to saying that the integral of a piece- 
wise discontinuous function (the integrand in Eq. I19|l is 
continuous. 



contrast to pi(x,t) and ph(x,t), however, the decay con- 
stant is not a monotonic function of time. This can be 
understood as follows: at early times, the large island- 
to-island distances come from large holes and therefore 
Pi2i{x,t) ~ ph{x), as mentioned earlier. (The inset of 
Fig. IBfc) confirms this.) However, as the island frac- 
tion f(t) approaches unity, the system becomes mainly 
covered by large islands, and Pi2i{x,t) should approach 
~ Pi(x,t) 2 asymptotically (see the second term in the 
bottom Eq.H| . 

In Fig. 8, we plot the decay constants for the three 
different distributions, th, t,, and r^,;. Note that when 
/ < 0.5, Th ~ U2i, as discussed above. As / — > 1, the 
behavior of Tin is controlled by ti , as suggested by Eq. 
20. Because pi2i ~ p\, we expect Tin — ► 0.5^; however, 



V. CONCLUSION 



To summarize, we have extended the KJMA model 
to the case where the homogeneous nucleation rate is 
an arbitrary function I(t) of time, deriving a number 
of analytic results concerning the properties of various 
domain distributions. We have also presented highly 
efficient simulation algorithms for ID nucleation-growth 
problems. Both analytical and simulation results are 
in excellent agreement. In the companion paper, we 
discuss the application of these results to experiments in 
general and to the analysis of DNA replication kinetics 
in particular. 
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